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EQUIVALENT SOURCE MODELING OF 
THE WAIN FIELD USING MAGSAT DATA 

M. A. Mayhew and R. H. Estes 
BTS, Inc. Seabrook, MO 20706 


Abstract 


An iterative least squares estimation algorithm with the capability 
for including a priori statistical information has been implemented to 
recover multiple magnetic dipole models of the Earth’s main magnetic 
field. The dipoles are fixed to a specified radius at or below the core- 
mantle boundary and centered on equal area blocks. The algorithm can 
solve for dipole magnitudes only (fixed orientations), or allow full 
freedom of orientation und solve for vector coniponents. External field 
parameters and observatory anomaly biases can also be estimated simultan- 
eously. Time dependence is modeled using first time derivatives for 
dipole vector components. Single-epoch and time dependent dipole models 
are derived using MAGSAT and observatory annual means data. Equivalent 
spherical harmonic representations are computed in closed fonn from the 
dipole models and compared with truncated spherical harmonic models 
estimated In the standard way from the same data sets. In particular, a 
21“ spatial resolution model based on 93 dipoles is computed based on 
observatory annual means data from 1974 through 1977 and a selected 
MAGSAT data set and is compared with candidate I6RF 1975 models and their 
1980 secular variation. The equivalent dipole source representation is 
shown to be comparable to the standard spherical harmonic approach in 
accuracy and, for high resolution models, to be superior in computational 
efficiency for field model evaluation when three deg ree-of- freedom uncon- 
strained dipoles are utilized. Fixing of the dipole positions results in 
rapid convergence of the dipole solutions for single-epoch models. For 
time dependent models, a sufficiently long time interval of data must be 
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available or a priori values and statistics available for the derivatives 
to achieve convergence. In contrast to standard spherical harmonic 
models based on the same data, the correlation structure of dipole 
magnetic moment derivatives with t'le constant magnetic moment parameters 
is small, indicating the ability t,o strongly separate the constant field 
and the secular variation. The dipole equivalent source approach for the 
main field was found not to be effective as a technique for 
distinguishing core and crustal fields. 
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1.0 INTRODUCTION 


The equivalent source technique is a conceptually simple means of 
describing potential fields which has been widely applied in exploration 
geophysics (e.g. Dampney, 1969), in modeling long wavelength satellite 
magnetic anomalies (Mayhew, 1979), and in modeling the core field (see 
review by Sterns and All dredge, 1973). The method consists of setting up 
an internal arrangement of dipoles and adjusting their magnetic moments 
(and in some applications their positions) such that the synthetic field 
arising from them collectively best fits a set of observations of the 
field. A number of authors have attempted to apply the equivalent source 
technique to modeling of the main field and its secular change. McNish 
(1940) and Lowes and Runcorn (1951) used graphical approximation 
techniques for this purpose. All dredge and Hurwitz (1964) and Alledredge 
and Sterns (1969) treated the problem using inverse theory and with the 
benefit of computers. The goal of the modern studies has been to gain 
insight into the physics of core field generation by examination of the 
behavior of modeled sources: for example, the fields of radial dipoles 

can be taken to approximate the fields of current loops. Some authors 
(e.g. Peddie, 1979) have used expressions for current loops directly. 

In most of these studies, the angular coordinates and radial 
distance of the dipoles are treated as variables in addition to the 
dipole moments, and herein lies the difficulty of this approach to core 
field modeling. The function expressing the dipole field is non-linear 
in the geographic variables, and this leads to a non-linear inverse 
problem. Many iterations through a given data set are required for 
convergence; for example, 180 iterations were required to fit 21 dipoles 
to the 1955 Finch-Leaton model of the field (Stearns and Al l dredge, 

1973), and 25,000 iterations were required to fit the 1966 IGRF secular 
change field (All dredge and Stearns, 1969). Furthermore, convergence was 
not achieved at all unless good starting parameter estimates were 
available. Previous approaches have sought to reproduce existing field 
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models expressed in spherical harmonics. A current function (Chapman and 
Bartels, 1940) is derived for a thin shell within the core from the 
spherical harmonic coefficients, and the starting positions of the 
dipoles are taken to be the current foci. The Inapplicability of this 
approach when one is trying to generate a field model from scratch Is 
obvious. Further, there Is no way of knowing before hand how many 
dipoles are required to properly model the field. 

Our previous work In magnetic dipole analysis has been aimed at 
modeling the long wavelength anomalies In the total magnetic field from 
satellite data* An array of dipoles Is laid out in equal area at the 
Earth's surface in a specific region, and the magnetic moments of the 
dipoles adjusted by a least squares procedure so as to collectively give 
rise to a field which best fits j;hat observed over the region. For this 
(essentially geologic) application, the dipoles are constrained to lie In 
the direction of the main field at the dipole position. In this study 
we have applied the methodology and the associated software to modeling 
the core field. We have experimented with fitting the main field using 
different numbers of equivalent sources at fixed radii at and within the 
core-mantle boundary. By fixing the radius for a given series of runs, 
we avoid the convergence problems which result from the extreme non- 
linearity of the problem when dipole positions are allowed to vary. The 
study had three main goals: 1) to compare the equivalent source approach 

with the standard spherical harmonic approach for modeling the main field 
In terms of accuracy and computational efficiency, 2) to examine the 
possibility of distinguishing core and crustal fields with this approach, 
and 3) to see whether a systematic distribution of solution parameters, 
possibly indicative of fluid motions, could be described, and if so what 
dipole density is required to resolve the pattern. 


\i 
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2.0 METHODOLOGY 

Mayhew (1978) and Mayhew and Estes (1980) have given the expressions 
for the anomaly components and the anomaly in the total field due to a 
dipole having arbitrary position and orientation at an arbitrary external 
position. The magnetic potential at an external point of a dipole of 

moment M is V = -M • V(l/A) , where A is the distance between the 

dipole and the external point; the associated vector field is F = - W. 
Let 0 be colatitude, <\> be longitude, and r be the radial distance, 
and let M have components (M^, M 0 , M^) . Then V may be written 

V = CMj,(rA-r') - MgrB + M^rC]/A^ , ( 1 ) 


where 


A = (r^ + r'^ - 

A = COS0 COS0* + sin0 sin0' cos(<|)-4'') 

B = COS0 sin0' + sin0 cos0' cos(<t>-<|i' ) 

C = sin0 sin(<(>-<|)' ) . 

Primed variables refer to the dipole position, unprimed to the external 
position. Then 
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The components of T are obtained by straightforward differentiation, 
and are complicated functions of the geographic coordinates and the 
components of the magnetic moment. If the dipole has a radial 
orientation the expressions for the components simplify; they are given 
by All dredge and Hurwitz (1964). 
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The goal of the analysis is to estimate values for the magnetic 
moments (M^, Mq, M^) for all dipoles from observations of the total 

field and its components. A software system (see Appendix) which is an 
elaboration of our previous equivalent source modeling software has been 
developed for this purpose. It uses an iterative least squares algorithm 
with the capability for a priori information which has been described in 
some detail by Mayhew and Estes (1980). The dipoles are set by user input 
to a specified radius within the earth, and are centered on nearly equal 
area blocks. The system has the option to solve for either the dipole 
magnitudes only (with the dipole orientations constrained to a preselected 
direction), or to allow full freedom* of orientation for the dipoles and 
solve for the three components. The procedure for obtaining an approximate 
equal area distribution of dipole locations was obtained by use of an 
icosahedron, a regular polyhedron with 20 equilateral triangular faces, 
inscribed in a sphere, with the edges radially projected onto the sphere 
forming “spherical polyhedrons." The spherical Icosahedron is the division 
of the sphere having the greatest number of regular pieces, and forms the 
base on which a nearly uniform distribution of points on the sphere may be 
defined. The following sets of points defining dipole locations have been 
implemented: 


number of points 12 42 92 162 

angular separation 64“ 32“ 21“ 16“ 
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When the dipole directions are unconstrained there are three parameters 
to be estimated per dipole, but only one per dipole when the direction is 
constrained radially, since in that case Mg = ■ 0. Secular variation 

is modeled by making each dipole magnetic moment component a linear 
function of time. This doubles the number of parameters for a given 
problem at a single epoch. An iterative technique is required for this 
problem when scalar magnetometer data is utilized due to the non- 
linearity introduced for this data type. The inversion problem is 
linear, however, for B^, Bg, B^ data. 

The set of simultaneously estimated parameters has been further 
extended in two ways. First, on option a set of three external field 
parameters expressed as spherical harmonic coefficients g’jQ, can 

be computed. Second, also on option, observatory anomaly biases can be 
estimated at each observatory when annual means data is used (Langel , 
Estes and Mead, 1982). 

The spherical harmonic representation of a given dipole distribution 
may be expressed in closed form. Let 


V(r,e,X) 




+ Pii(cos6)) 


( 3 ) 


be the scalar potential of the main field where the Legendre polynomials 
are Schmidt normalized. The last term represents that portion of the 
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potential originating from sources outside of the sphere of radius 11 
(three parameter external field model). The magnetic field is then 

T » -W . 

r , 

Consider the set of |j} dipoles, where the i dipole is described by a 
source vector 


and a spatial position 


Mr . Mfl , M. 

r, 9, 


r,, 9,, . 


The spherical harmonic coefficients are given (in Schmidt normalized 
form) by 

, -^(n+2)(j} n-1 , , 

ill 'i " "'•i 




sin m^ 


■ 0 . 30 

1 


cos mX - m P„(cos9,) 


(4) 


^ , .-(n+2) {^} n-1 , , 

''t n PJws 9,) sin mx. 


+ M 


9Pnm(c^s«i) 


0 ^ 30 


si 


cos mX 


n six, t n, P„„(cos9,) - 5 ^ 


The coefficients may be calculated up to any order. The spatial power 
spectral content of the expansion is given by the relation (Lowes, 1966, 
1974) 
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3.0 RESULTS 

Several main field solutions have been derived for various numbers 
of dipoles arranged on equal area projections at the ccre-mantle 
boundary. In these various solutions, the dipoles have been either fixed 
radially or allowed complete freedom of orientation, and have been 
modeled with and without time dependence and with and without the 
addition of a geocentric dipole (giving three additional parameters). 

3.1 Single-epoch Models 

Based on the MGST(6/80) data set of quiet November 5 and 6, 1979, 
MAGSAT data (Langel et a1 , 1980), several dipole models Identified In 
Table I have been derived. The models for which the dipoles are 
constrained radially have a single degree of freedom’ per dipole, while 
the unconstrained models have three degrees of freedom per dipole. The 
inclusion of the geocentric dipole significantly Improved the rms error 
of the 32° dipole density solution. .The direction of the geocentric 
dipole moment from the model was within 10° of that calculated from the 
first degree spherical harmonic terms of MGST(6/80), while the magnitude 
was greater by approximately The magnitude of the solution 
geocentric dipole dominated those of the core/mantle boundary by an order 
of magnitude. The Inclusion of the geocentric dipole In the 21° dlpolJ 
density solution of model 4 showed no improvement in the fit to the data 
over model 3 and the geocentric dipole showed no closer agreement to the 
first degree terms of MGST(6/80) than did that of the 32° dipole density 
solution. All parameters In these solutions were strongly observable In 
the least squares estimation, with standard errors from three to four 
orders of magnitude smaller than the magnitude of the parameters. The 
coefficients of the spherical harmonic representation given by Equation 
(4) for models 3, 4 and 5 from Table I showed differences of less than 1 
nT (nanotesla) when compared with MGST(6/80) though degree and order 
ten. The spatial power spectra computed from spherical harmonic 
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expansions of the d1po(e models were analyzed relative to crustal and 
core content. Figure I displays the power spectra for one of the dipole 
models and M6ST(10/81), a spherical harmonic model based on MAGSAT data 
through degree and order 23 (Langel and Estes, 1982a). 

The cause of the apparent bulge In the power spectrum for Dipole 
Model 14 (Figure 1) was Investigated via a simulation. The simulation 
consisted of generating synthetic data at exactly the same spatial 
locations as the MGST(6/80) data set used 1n the equivalent dipole source 
models (a "pseudo-data set"), with the "data" consisting of measurements 
due to a crustal field plus a geocentric dipole only. The crustal field 
model used was a global equivalent source model based on POGO data 
(R. A. Langel, personal communication) which was expanded to degree and 
order 40 In spherical harmonics to analyze Its spectral character, as 
well as to generate the synthetic data. The spectrum of the synthetic 
crustal field Is undoubtedly biased at low order, since ultra-long wave- 
length trends were removed from the POeO data In generating the crustal 
equivalent source model. The work of Meyer et al (pvvsonal communica- 
tion) predicts a white spectrum down to very lofe order. The distortion 
of the equivalent source spectrum does not, however. Invalidate certain 
Inferences which can be made from the results of the simulation. Least 
squares fits were made to the pseudo-data set with a single geocentric 
dipole and with a geocentric dipole plus dipoles of 32®x32° and 21*x21‘’ 
density at the core/mantle boundary. The results for cares with and 
without noise are given In Table II, and the spectra of selected results 
are given in Figure 2. The results seem to Indicate that the beginning 
of the bump In the spectrum of Dipole Model #4 In Figure 1 Is due to 
crustal Influence, while the departure of the spectrum from that of 
MGST(10/81) around expansion order 17 Is due to the resolution limits of 

the dipole density (I.e., =17). 
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Time-dependent Hodels 

Dipole solution models have been obtained with both 32° and 21° 
resolutions based on a quiet (Kp<l") MAGSAT data set extending over four 
months (November 5, 1979 - March 15, 1980). Time dependence was modeled 
using first time derviatives for the dipole magnetization vector 
components, which doubles the total number of parameters in the 
solution. The solutions displayed va very slow convergence in the time 
derivatives (the derivatives had very large magnitudes), although at each 
iteration the conversion of the dipole parameters to spherical harmonic 
coefficients g^,„ , h^^^ and g^^j^ , by equation (4) showed close agree- 
ment with a standard spherical harmonic model to degree and order 13 in 
the constant terms and degree and order 7 in the first derivative terms 
based on the same data set. The observability of the time variation of 
the individual dipoles over the short time period was poor, while the 
information content to the extent of determining the global time 
variation to degree 7 in spherical harmonics was strong. A priori values 
for the geocentric dipole derivatives were obtained from g^^Q , g^j^ and 
hj^j^ of the standard spherical harmonic model and the solution was then 
statistically constrained for the geocentric dipole using these values. 
The convergence using this technique showed much Improvement. 


ti 

1 








To improve the observability of the time variation, observatory 
annual means data for a selected set of magnetic observatories from 
1974-1977 was added to the above MAGSAT data set, and an unconstrained 
epoch 1980 dipole solution of 21° resolution (with a geocentric dipole) 
was obtained splifing simultaneously for the observatory anomaly vectors 
(Langel et al , 1982). This solution is referred to as the 1974-1980 
Dipole Model. The solution converged in three iterations, showing small 
correlations among the constant and derivative parameters (typically on 
the order of .01), although a few correlations reached values of 
approximately .4 between some nearest neighbor dipole parameters. The 
derivatives were on the order of one to two orders of magnitude smaller 
than the constant dipole parameters, with standard errors on the same 
order as the derivatives. This correlation structure is in stark 
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contrast to that observed in standard spherical harmoKic models based on 
the same data set where the constant and derivative terms for a given 
degree and order may show correlations as high as .99. The dipole 
magnetization values per unit volume and their derivatives are given in 
Table III. With the exception of the geocentric dipole, all dipoles were 
placed at a geocentric radius of 3000 km. The volume associated with 
each dipole was that of a spherical shell of 40 km, thickness. 

Volume « (21® X x 3000 km)^ x 40 km . 


The same volume was arbitrarily assigned to the geocentric dipole. The 
observatory anomaly bias vectors are presented in Table IV, while the 
coefficients of the spherical harmonic expansion of the 93 dipoles 
through degree and order 24 (as calculated from equation (4) using mag- 
netic moments given by dipole magnetizations multiplied by their 
respective volumes) are given in Table V. 

The 1974-1980 Dipole Model is compared with candidate models for the 

1975 IGRF presented by NASA, USGS and IGS in Tables VI through VIII. 

Table VI evaluates the fit to selected observatory annual means (Langel 

and Estes, 1982b) using the statistic a defined by Langel et al . (1982), 

while Table VII compares the models using MAGSAT data. The dipole model 

performs well with respect to the candidate IGRF models. Table VIII 

displays the low degree and order spherical harmonic coefficients of the 

Dipole Model and the USGS candidate IGRF Secular Variation model, which 

is based on conventional techniques using more recent observatory annual 

means. The trends in the coefficients are very similar, although the 

magnitudes of the first degree dipole model terms are larger by a few nT 

than the USGS model. This is a result of the dipole model heavily 

utilizing Magsat data at epoch 1980. Field models based only on Magsat 
• • 
data yield gio on the order of 28 nT/yr and hii on the order of 

-20 nT/yr. The availability of ample observatory annual means data for 

years later than 1980 are required to more fully evaluate the different 

secuiar variation terms. 


10 







^ j 







Business. iNi) TEaimimicii. Systems, !m 


4.0 CONCLUSIONS 

The modeling of the main field with an equivalent dipole 
representation has been found to be comparable to the standard spherical 
harmonic approach in accuracy. The 32® dipole density (42 dipoles) 
corresponds approximately to an eleventh degree/order spherical harmonic 
expansion (143 parameters) while the 21® dipole density (92 dipoles) 
corresponds to approximately a seventeenth degree and order expansion 
(288 parameters). Comparison of the number of arithmetic operations 
(multiplications, additions and trigonometric and square root 
evaluations) required to evaluate the magnetic field from dipole models 
and spherical harmonic models shows the dipole models to be superior in 
computational efficiency for spatial resolution better than 21® whon 
three degrees-of- freedom unconstrained dipoles are utilized. The 
computational burden of field model evaluation is comparable to the 
standard spherical harmonic approach when the dipoles are constrained 
radially. Fixing of the dipole positions results in rapid convergence of 
the dipole solutions fo'* single-epoch models. For time dependent models, 
a sufficiently long time interval of data or a priori values and 
statistics must be available for the derivatives to achieve convergence. 
In contrast to standard spherical harmonic models based on the same data, 
the correlation structure of dipole magnetic moment derivatives with the 
constant magnetic moment parameters is small, indicating the ability to 
strongly separate the constj^nt field and the secular variation. The 
dipole model of 21° spatial resolution using observatory annual means 
data from 1974-1977 together with MAGSAT data performs very well when 
compared to candidate IGRF models of the main field and its secular 
variation based on conventional spherical harmonic techniques. 

The dipole equivalent source approach for the main field was found 
not to be effective as a technique for distinquishing core and crustal 
fields. The investigation of possible indications of core fluid motions 
available from the equivalent source dipole model is left for further 
study. 
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FIGURE 1 Power spectrum of field model MGST (10/81) 
compared with spectra of Dipole Models 4 and 5 
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FIGURE 2 Power spectra of synthetic "crustal" field and 

two dipole array models fit to this field compared with 
spectra of field model MGST (10/81) and Dipole Model 4 ■; 
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TABLE III 


DIPOLE MAGNETIZATION VALUES PER UNIT VOLUME 
FOR 1974-1980 DIPOLE MODEL 
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TABLE VII STATISTICAL COMPARISON OF CANDIDATE 
IGRF MODELS AND 1974-1980 DIPOLE MODEL 
ON SELECTED MAGSAT DATA (nT) AT EPOCH 1980 


Spherical 

Harmonic Expansion 



Order 

"b 

"x 

"y 

"z 

NASA 

(10) 

13.3 

11.3 

11.2 

14.9 

USGS 

(10) 

15.7 

10.3 

10.2 

13.0 

IGS 

(10) 

112.0 

85.2 

66.9 

118.0 

DIPOLE 

(16) 

7.0 

6.8 

8.0 

8.6 

DIPOLE 

(10) 

11.2 

10.0 

10.2 

12.4 
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Table VIII 

Low degree and order spherical siarmonlc coefficients at Epoch 1980 
derived from EQUIVALENT DIPOLE MODEL of 21® resolution using MAGSAT data 
and observatory Annual Means from 1974-1977. 


USGS IGRF 
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Dipole Model Model Epoch 1980 


Degree 

Order 

®nm 

^nm 

^nm 

linm 

nm 

®nm 
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APPENDIX 

SOFfWARE DESCRIPTION AND FLOW CHART 


The software system was developed to recover a multiple dipole 
representation of the main field and is coded in FORTRAN. An assembly 
language subroutine is utilized to obtain required computer core memory 
dynamically in the IBM 360 and 3081 environment. The analysis utilizes 
an iterative least squares algorithm with the capability for a priori 
information. The dipoles are set by user input to a specified radius 
within the earth, and are centered on equal area blocks. The system has 
the option to solve for either the dipole magnitudes only (with the 
dipole orientations constrained to a preselected direction), or to allow 
full freedom of orientation for the dipoles and solve for the three 
components . 

The least squares equations are as follows: 


6Xn+i = (A^WA + 


A^W6y„ 




where 

A is the partial derivative matrix of the measurements 
A respect to the parameters 

X is the vector of adjusted dipole parameters 

6y is the vector of residuals 

W is the weight matrix for the measurements 

Aq is the a priori parameter covariance matrix 

A 

Xq is the a priori estimate of the parameters and the 
estimate at the (n+1)^^ iteration is: 

Vl ' • 


A-1 
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An iterative technique is required for this problem when scalar 
magnetometer data is utilized due to the non-linearity introduced for this 
data type. The inversion problem is linear, however, for B^, Bg, B^ data. 


The major program subroutines and their functions are as follows: 

MAIN: Program driver. Calls routines to initialize all variables. 

Loops through the least squares iteration, calling for data 
and parti als, and accumulates the normal matrix and right-hand- 
side. Calls for inversion of normal matrix. Calculates updated 
parameter vector and statistics. After last iteration calculates 
statistics of data set for the recovered solution. 

TSINV: Inverts a symmetric positive definite matrix stored row-wise in 

upper triangular form. 

A A 

SETRHS: Initializes the right-hand-side to A"^(Xq - for each itera- 

tion when a priori information is used. 

CLEAR: Sets arrays to zero values. 

FUN: Canputes the partial derivatives of B^,, Bg, B^ and B with respect 

to the parameters and calculates values for B^, Bg, B^, B in 
terms of current estimates of the parameters. 

FUN2: Entry in FUN. Initializes the positions of dipoles on first 

call, establishing the size of equal area blocks to be used. 
Zeroes the normal matrix and right-hand-side on each call. 
Calculates the a priori covariance matrix when there are 
constraints for the orientation of the dipoles along specified 
di recti ons. 


A-2 
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FUN3: Entry in FUN. Calculates the angles between the solution 

dipoles and specified directions. 

TMULT: Calculates the matrix B x B^. 



FLD; Specifies a preferred direction for the dipoles. 



!■[ 


[1 


; : > ■ I 

r; ] 

'ii 



EQAREA: Positions the dipoles at the center of equal area blocks with 

optional angular separations of 64®, 32®, 21®, and 16®. 

GNORML: Computes random numbers of specified mean and variance. 

CORLPR: Prints solution normal, covariance and correlation matrices. 

SVDATA: Inputs MAGSAT scalar and vector data from tape. Provides data 

®r’ ^0’ ^ 41 ’ ® ^ particular time and position on each call. 

SPHCOE: Computes spherical harmonic coefficients from the dipole 

parameters. 

GSPACE: Assembly language routine which computes machine core memory 

size requirements dynamically. 

BIAS: Computes the observatory anomaly biases on option w'nen obser- 

vatory annual means data is used. 

SDATA: Inputs observatory and repeat data from tape. 


I !T'* 
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Specify prpgram options 
in data statements 


Initialize 

statistical I 




SET Number of least 
squares Iterations, ITER 


IT - IT + 1 


CALL FUN 2 


CALL SVOATA 


ALL FUN 


Accumulate statistical 
variables 


lt4»ll’W^TOPrrFlt 


form partial derivatives 


li Accumulate normal matrix | 

P and right-hand-side ! 

r ii 

(^ST DATA P0IN)--{^ 

^ i 


! m& 
.] 1 

CALL TSINV 
— ■ — 


WIHTmNiMnTOJI 


complete iteration 
statistics and display 


LAST ITERATION 


MM.i *in>f 


through data set and displa 



















